#ifndef SPHERE_H
#define SPHERE_H

#include "./sMathLib/Vector3/Vector3.hpp"
#include "Box3D.h"
#include <cmath>
#include <map>
#include <vector>
#include <iostream>

class NJRvector3d;

namespace ibm
{

template<typename T, template<typename U> class BaseObject>
class Sphere: public BaseObject<T>
{
protected:

	T volume;

	bool volumeCached;

	Box3D<T> boundingBox;
	
	std::vector< sMathLib::VectorX<T, 3> > shiftedPositions;

	std::vector< Box3D<T> > shiftedBoundingBox;
	

public:

	Sphere();

	Sphere(sMathLib::VectorX<T, 3>& inputPosition, sMathLib::VectorX<T, 3>& inputVelocity, sMathLib::VectorX<T, 3>& inputAVelocity, T& inputRadius);

	Sphere(NJRvector3d& inputPosition, NJRvector3d& inputVelocity, NJRvector3d& inputAVelocity, T& inputRadius);

	Sphere(T& px, T& py, T& pz, T& ux, T& uy, T& uz, T& avx, T& avy, T& avz, T& inputCharaLength);

	Sphere(const Sphere<T, BaseObject>& rhs);

public:

	const T& getRadius() const;

	const T& getVolume() const;
		
	void createBoundingBox(T& driftIn, T& driftOut);
	
	void createBoundingBoxForShiftedSphere(T& driftIn, T& driftOut);
		
	void generateShiftedSpheres(Box3D<T>& computationDomain);
	
	unsigned getShiftedSize();
	
	sMathLib::VectorX<T, 3>& getShiftedPosition(unsigned index);
	
	Box3D<T>& getShiftedBounding(unsigned index);



		

};

template<typename T, template<typename U> class BaseObject>
Sphere<T, BaseObject>::Sphere()
	:volumeCached(false)
{

}

template<typename T, template<typename U> class BaseObject>
Sphere<T, BaseObject>::Sphere(sMathLib::VectorX<T, 3>& inputPosition, sMathLib::VectorX<T, 3>& inputVelocity, sMathLib::VectorX<T, 3>& inputAVelocity, T& inputRadius)
	:BaseObject(inputPosition, inputVelocity, inputAVelocity, inputRadius), volumeCached(false), volumeCached(false)
{

}

template<typename T, template<typename U> class BaseObject>
Sphere<T, BaseObject>::Sphere(NJRvector3d& inputPosition, NJRvector3d& inputVelocity, NJRvector3d& inputAVelocity, T& inputRadius)
	:BaseObject<T>(inputPosition, inputVelocity, inputAVelocity, inputRadius), volumeCached(false)
{

}

template<typename T, template<typename U> class BaseObject>
Sphere<T, BaseObject>::Sphere(T& px, T& py, T& pz, T& ux, T& uy, T& uz, T& avx, T& avy, T& avz, T& inputCharaLength)
	:BaseObject<T>(px, py, pz, ux, uy, uz, avx, avy, avz, inputCharaLength), volumeCached(false)
{

}

template<typename T, template<typename U> class BaseObject>
Sphere<T, BaseObject>::Sphere(const Sphere<T, BaseObject>& rhs)
	:BaseObject<T>(rhs)
{
	*this = rhs;
}

template<typename T, template<typename U> class BaseObject>
const T& Sphere<T, BaseObject>::getRadius() const
{
	return BaseObject<T>::getCharacteristicLength();
}

template<typename T, template<typename U> class BaseObject>
const T& Sphere<T, BaseObject>::getVolume() const
{	 
	if(!volumeCached)
		volume = (4. / 3.) * (4 * std::atan(1.)) * (getRadius() * getRadius() * getRadius());
	return volume;
}

template<typename T, template<typename U> class BaseObject>
void Sphere<T, BaseObject>::generateShiftedSpheres(Box3D<T>& computationalDomain)
{
	shiftedPositions.push_back(this->position);
	
	typedef sMathLib::VectorX<T, 3> Vector3;	

	if(boundingBox.x0 < computationalDomain.x0)
	{		
		int Lx = computationalDomain.x1 - computationalDomain.x0;
		shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y(), BaseObject<T>::getPosition().z() ) );
	
		if(boundingBox.y0 < computationalDomain.y0)
		{
			int Ly = computationalDomain.y1 - computationalDomain.y0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z() ) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z() ) );
				
			if(boundingBox.z0 < computationalDomain.z0)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()+Lz) );
			}
			else if(boundingBox.z1 > computationalDomain.z1)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()-Lz) );
			}
			
		}
		else if(boundingBox.y1 > computationalDomain.y1)
		{
			int Ly = computationalDomain.y1 - computationalDomain.y0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z() ) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z() ) );
			if(boundingBox.z0 < computationalDomain.z0)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()+Lz) );
			}
			else if(boundingBox.z1 > computationalDomain.z1)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()+Lx, BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()-Lz) );
			}
			
		}
	
	}
	else if(boundingBox.x1 > computationalDomain.x1)
	{
		int Lx = computationalDomain.x1 - computationalDomain.x0;
		shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y(), BaseObject<T>::getPosition().z() ) );
	
		if(boundingBox.y0 < computationalDomain.y0)
		{
			int Ly = computationalDomain.y1 - computationalDomain.y0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z() ) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z() ) );
			
			if(boundingBox.z0 < computationalDomain.z0)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()+Lz) );
			}
			else if(boundingBox.z1 > computationalDomain.z1)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()-Lz) );
			}
			
		}
		else if(boundingBox.y1 > computationalDomain.y1)
		{
			int Ly = computationalDomain.y1 - computationalDomain.y0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z() ) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z() ) );
			
			if(boundingBox.z0 < computationalDomain.z0)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()+Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()+Lz) );
			}
			else if(boundingBox.z1 > computationalDomain.z1)
			{
				int Lz = computationalDomain.z1 - computationalDomain.z0;
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()   , BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()-Lz) );
				shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x()-Lx, BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()-Lz) );
			}
			
		}
		
	}	

	
	else if(boundingBox.y0 < computationalDomain.y0)
	{
		int Ly = computationalDomain.y1 - computationalDomain.y0;
		shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z() ) );
		
		if(boundingBox.z0 < computationalDomain.z0)
		{
			int Lz = computationalDomain.z1 - computationalDomain.z0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()+Lz) );
		}
		else if(boundingBox.z1 > computationalDomain.z1)
		{
			int Lz = computationalDomain.z1 - computationalDomain.z0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()+Ly, BaseObject<T>::getPosition().z()-Lz) );
		}
	}
	else if(boundingBox.y1 > computationalDomain.y1)
	{
		int Ly = computationalDomain.y1 - computationalDomain.y0;
		shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z() ) );
		
		if(boundingBox.z0 < computationalDomain.z0)
		{
			int Lz = computationalDomain.z1 - computationalDomain.z0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()+Lz) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()+Lz) );
		}
		else if(boundingBox.z1 > computationalDomain.z1)
		{
			int Lz = computationalDomain.z1 - computationalDomain.z0;
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
			shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()-Ly, BaseObject<T>::getPosition().z()-Lz) );
		}
	}
	else if(boundingBox.z0 < computationalDomain.z0)
	{
		int Lz = computationalDomain.z1 - computationalDomain.z0;
		shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y(), BaseObject<T>::getPosition().z()+Lz) );
	}
	else if(boundingBox.z1 > computationalDomain.z1)
	{
		int Lz = computationalDomain.z1 - computationalDomain.z0;
		shiftedPositions.push_back( Vector3( BaseObject<T>::getPosition().x(), BaseObject<T>::getPosition().y()   , BaseObject<T>::getPosition().z()-Lz) );
	}
}

template<typename T, template<typename U> class BaseObject>
void Sphere<T, BaseObject>::createBoundingBoxForShiftedSphere(T& driftIn, T& driftOut)
{
	int r_extend = (int)( getRadius() + std::abs(driftOut) ) + 1;
	
	int d_extend = 2 * r_extend;
	
	for(unsigned sphereIndex = 0; sphereIndex != shiftedPositions.size(); ++sphereIndex)
	{
		shiftedBoundingBox.push_back( Box3D<T>() );
		Box3D<T>& currentBounding = shiftedBoundingBox[sphereIndex];
		const sMathLib::VectorX<T, 3>& currentPosition = shiftedPositions[sphereIndex];
			
		currentBounding.x0 = static_cast<int>(currentPosition[0]) - r_extend;
		currentBounding.x1 = currentBounding.x0 + d_extend;
			
		currentBounding.y0 = static_cast<int>(currentPosition[1]) - r_extend;
		currentBounding.y1 = currentBounding.y0 + d_extend;
			
		currentBounding.z0 = static_cast<int>(currentPosition[2]) - r_extend;
		currentBounding.z1 = currentBounding.z0 + d_extend;
	}
}

template<typename T, template<typename U> class BaseObject>
void Sphere<T, BaseObject>::createBoundingBox(T& driftIn, T& driftOut)
{
	int r_extend = (int)( getRadius() + std::abs(driftIn) ) + 1;
	
	int d_extend = 2 * r_extend;
	
	boundingBox.x0 = static_cast<int>(BaseObject<T>::getPosition()[0]) - r_extend;
	boundingBox.x1 = boundingBox.x0 + d_extend;
	
	boundingBox.y0 = static_cast<int>(BaseObject<T>::getPosition()[1]) - r_extend;
	boundingBox.y1 = boundingBox.y0 + d_extend;
	
	boundingBox.z0 = static_cast<int>(BaseObject<T>::getPosition()[2]) - r_extend;
	boundingBox.z1 = boundingBox.z0 + d_extend;
	
}

template<typename T, template<typename U> class BaseObject>
sMathLib::VectorX<T, 3>& Sphere<T, BaseObject>::getShiftedPosition(unsigned index)
{
	return this->shiftedPositions[index];
}

template<typename T, template<typename U> class BaseObject>
unsigned Sphere<T, BaseObject>::getShiftedSize()
{
	return shiftedPositions.size();
}

template<typename T, template<typename U> class BaseObject>
Box3D<T>& Sphere<T, BaseObject>::getShiftedBounding(unsigned index)
{
	return this->shiftedBoundingBox[index];
}

//

} /* namespace ibm */
#endif /* SPHERE_H */
